function [out, vars, vectors] = GPOPSRobustnessComputeVars(output,setup)
% computes variables based on GPOPS result, including various
% statistics, and terms from the optimal tax formula. Also does
% some plotting, but that is not the main purpose. Stores
% transformed results.
    
tic
%% assign variables
    auxdata = setup.auxdata;
    
    %% Parameters
    inpscale = auxdata.settings.inpscale;
    out.L_M = inpscale.L_M.* output.result.solution.parameter(1,1);
    out.L_R = inpscale.L_R.* output.result.solution.parameter(1,2);
    out.L_C = inpscale.L_C.* output.result.solution.parameter(1,3);
    out.K_B_M = inpscale.K_B_M.* output.result.solution.parameter(1,4);
    out.K_B_R = inpscale.K_B_R.* output.result.solution.parameter(1,5);
    out.K_B_C = inpscale.K_B_C.* output.result.solution.parameter(1,6);
    out.K_E = inpscale.K_E.* output.result.solution.parameter(1,7);
    out.K_S = inpscale.K_S.* output.result.solution.parameter(1,8);
    out.b = inpscale.transfer.* output.result.solution.parameter(1,9);

    out.K_B = out.K_B_M + out.K_B_R + out.K_B_C;
    
    % derived factor prices
    gpops_fun = auxdata.gpops_fun;

    out.Y = gpops_fun.prod_fun.Y(out);
    out.Y_M = gpops_fun.prod_fun.Y_M(out);
    out.Y_R = gpops_fun.prod_fun.Y_R(out);
    out.Y_C = gpops_fun.prod_fun.Y_C(out);
    out.Y_B_M = gpops_fun.prod_fun.Y_B_M(out);
    out.Y_B_R = gpops_fun.prod_fun.Y_B_R(out);
    out.Y_B_C = gpops_fun.prod_fun.Y_B_C(out);
    out.Y_E = gpops_fun.prod_fun.Y_E(out);
    out.Y_S = gpops_fun.prod_fun.Y_S(out); 

    
    out.wtilde = output.result.interpsolution.phase.time;
    out.wtilde_orig = out.wtilde;
    out.wtilde_noninterp = output.result.solution.phase.time;
    
    if auxdata.settings.log_w == 1
        if auxdata.settings.standardize == 1
            out.w_orig = exp(out.wtilde_orig.*auxdata.sigma_fit+auxdata.mu_fit);
            out.w = exp(out.wtilde.*auxdata.sigma_fit+auxdata.mu_fit);
            out.w_noninterp = exp(out.wtilde_noninterp.*auxdata.sigma_fit+auxdata.mu_fit);
        else
            out.w_orig = exp(out.wtilde_orig);
            out.w = exp(out.wtilde);
            out.w_noninterp = exp(out.wtilde_noninterp);
        end
    end
    
    out.w0 = auxdata.w0;
    out.wf = auxdata.wf;

    out.u = output.result.interpsolution.phase.state .* inpscale.u;
    out.u_noninterp = output.result.solution.phase.state .* inpscale.u;
    out.l = output.result.interpsolution.phase.control(:,1) .* inpscale.l;
    out.l_noninterp = output.result.solution.phase.control(:,1) .*inpscale.l;

    % multipliers
    out.mu = -output.result.solution.phase.integralmultipliers(5);
    out.xi_M = -output.result.solution.phase.integralmultipliers(6)./out.mu*out.Y_M;
    out.xi_R = -output.result.solution.phase.integralmultipliers(7)./out.mu*out.Y_R;
    out.xi_C = -output.result.solution.phase.integralmultipliers(8)./out.mu*out.Y_C;
    out.eta = output.result.interpsolution.phase.costate/out.mu;
    
    %% finer interpolation
    out.l_orig = out.l;
    out.u_orig = out.u;
    out.eta_orig = out.eta;
    
    out.l = interp1(out.w_orig,out.l_orig,out.w);
    out.u = interp1(out.w_orig,out.u_orig,out.w);
    out.eta = interp1(out.w_orig,out.eta_orig,out.w);
    
    %% assign parameters
    epsilon = gpops_fun.util_fun.epsilon;
    k = 1/epsilon;
    
    xi = gpops_fun.util_fun.xi;
    %% construct c
    out.c = 1./xi.*out.l.^(1+1/epsilon)/(1+1/epsilon) + out.u;
    
    %% compute derived variables
    out.mtax = 1-1./out.w.*1./xi.*out.l.^k;

    out.y = out.l.*out.w;
    out.tax = out.y - out.c;
    out.avg_tax = out.tax./out.y;
    
    % total mass
    out.total_mass = auxdata.gpops_fun.skill_dist.mass;
    
    % capital income and capital expenditure
    out.capital_income = out.K_B_M .* out.Y_B_M + ...
                out.K_B_R .* out.Y_B_R + ...
                out.K_B_C .* out.Y_B_C + ...
                out.K_E .* out.Y_E + ...
                out.K_S .* out.Y_S;

    out.q_B = gpops_fun.q_B;
    out.q_E = gpops_fun.q_E;
    out.q_S = gpops_fun.q_S;
    out.prod_fun_eta_R = auxdata.gpops_fun.prod_fun.eta_R;
    out.prod_fun_kappa_R = auxdata.gpops_fun.prod_fun.kappa_R;
    
    out.capital_expenditure = out.K_B_M .*out.q_B + ...
                out.K_B_R .* out.q_B + ...
                out.K_B_C .* out.q_B + ...
                out.K_E .* out.q_E + ...
                out.K_S .* out.q_S;
    
    % tax revenue and value of resource constraints
    out.labor_income = out.L_M * out.Y_M + ...
        out.L_R * out.Y_R + ...
        out.L_C * out.Y_C;

    % income shares
    out.labor_income_share = out.labor_income/out.Y;
    out.robot_income = out.K_B * out.Y_B_M; % requires equal
                                            % returns to all robot capital
    out.robot_income_share = out.robot_income/out.Y;
    out.labor_share_in_labor_plus_robot_income = out.labor_income/ ...
        (out.labor_income + out.robot_income);
    
    % bounds check
    fprintf('solution - lower bound, should be positive\n')
    check_lb = output.result.solution.parameter - setup.bounds.parameter.lower
    
    fprintf('upper bound - solution, should be positive\n')
    check_ub = setup.bounds.parameter.upper - output.result.solution.parameter
    
    % taxes
    out.tau_B = out.Y_B_M/out.q_B - 1;
    out.tau_E = out.Y_E/out.q_E - 1;
    out.tau_S = out.Y_S/out.q_S - 1;
    
    out.tau_return_B = 1 - 1/(1+out.tau_B);
    out.tau_return_E = 1 - 1/(1+out.tau_E);
    out.tau_return_S = 1 - 1/(1+out.tau_S);

    out.fval = output.result.objective * inpscale.welfare;
    out.welfare =  -output.result.objective * inpscale.welfare;
    
    out.welfare_M = output.result.solution.phase.integral(1);
    out.welfare_R = output.result.solution.phase.integral(2);
    out.welfare_C = output.result.solution.phase.integral(3);
    out.welfare_nonpart = output.result.solution.phase.integral(4);
    
    %% population shares, participation and employment
    w_lb = out.w(1);
    w_ub = out.w(end);
    
    % approximate phi_tilde
    out.phi_tilde_w = out.u - out.b;
    
    phi_tilde_fun = @(w) interp1(out.w,out.phi_tilde_w,w,'spline','extrap');
    
    w_range = w_lb:w_ub;
    figure(7)
    plot(out.w,out.phi_tilde_w,w_range,phi_tilde_fun(w_range),'--')
    title('phi tilde')
    legend('actual','approximated')
    
    
    f_M = @(w) gpops_fun.skill_dist.f_M(w,out.Y_M,out.Y_R,out.Y_C);
    f_R = @(w) gpops_fun.skill_dist.f_R(w,out.Y_M,out.Y_R,out.Y_C);
    f_C = @(w) gpops_fun.skill_dist.f_C(w,out.Y_M,out.Y_R,out.Y_C);

    out.f_M = f_M(out.w);
    out.f_R = f_R(out.w);
    out.f_C = f_C(out.w);
    out.f = out.f_M + out.f_R + out.f_C;
    
    out.pop_M = integral(f_M,w_lb,w_ub);
    out.pop_R = integral(f_R,w_lb,w_ub);
    out.pop_C = integral(f_C,w_lb,w_ub);
    
    if auxdata.settings.fix_participation
        
        pi_M = auxdata.settings.pi_M;
        pi_R = auxdata.settings.pi_R;
        pi_C = auxdata.settings.pi_C;
    
        h_M = @(w) gpops_fun.skill_dist.f_M(w,out.Y_M,out.Y_R,out.Y_C).*pi_M;
        h_R = @(w) gpops_fun.skill_dist.f_R(w,out.Y_M,out.Y_R,out.Y_C).*pi_R;
        h_C = @(w) gpops_fun.skill_dist.f_C(w,out.Y_M,out.Y_R,out.Y_C).*pi_C;
    
    else
        h_M = @(w) gpops_fun.skill_dist.f_M(w,out.Y_M,out.Y_R,out.Y_C).*gpops_fun.part_dist.G_M(phi_tilde_fun(w));
        h_R = @(w) gpops_fun.skill_dist.f_R(w,out.Y_M,out.Y_R,out.Y_C).*gpops_fun.part_dist.G_R(phi_tilde_fun(w));
        h_C = @(w) gpops_fun.skill_dist.f_C(w,out.Y_M,out.Y_R,out.Y_C).*gpops_fun.part_dist.G_C(phi_tilde_fun(w));

    end
    
    out.h_M = h_M(out.w);
    out.h_R = h_R(out.w);
    out.h_C = h_C(out.w);
    out.h = out.h_M + out.h_R + out.h_C;    

    out.F_M = zeros(size(out.w));
    out.F_R = zeros(size(out.w));
    out.F_C = zeros(size(out.w));
    out.H_M = zeros(size(out.w));
    out.H_R = zeros(size(out.w));
    out.H_C = zeros(size(out.w));

    % can use G_M here since participation costs are participation-independent
    out.pi_w = gpops_fun.part_dist.G_M(phi_tilde_fun(out.w)); 
    
    for i=1:length(out.H_M)
        out.F_M(i) = integral(f_M,w_lb,out.w(i));
        out.F_R(i) = integral(f_R,w_lb,out.w(i));
        out.F_C(i) = integral(f_C,w_lb,out.w(i));
        
        out.H_M(i) = integral(h_M,w_lb,out.w(i));
        out.H_R(i) = integral(h_R,w_lb,out.w(i));
        out.H_C(i) = integral(h_C,w_lb,out.w(i));
    end
    
    out.F = out.F_M + out.F_R + out.F_C;
    out.H = out.H_M + out.H_R + out.H_C;
    
    out.part_M = integral(h_M,w_lb,w_ub);
    out.part_R = integral(h_R,w_lb,w_ub);
    out.part_C = integral(h_C,w_lb,w_ub);
    
    out.part = out.part_M + out.part_R + out.part_C;
    
    % average participation rates by occ.
    out.pi_M = out.part_M/out.pop_M;
    out.pi_R = out.part_R/out.pop_R;
    out.pi_C = out.part_C/out.pop_C;
    
    out.h_M_cond_part = out.h_M / out.part_M;
    out.h_R_cond_part = out.h_R / out.part_R;
    out.h_C_cond_part = out.h_C / out.part_C;
    out.h_cond_part = out.h/out.part;

    out.H_M_cond_part = out.H_M / out.part_M;
    out.H_R_cond_part = out.H_R / out.part_R;
    out.H_C_cond_part = out.H_C / out.part_C;
    out.H_cond_part = out.H/out.part;
    
    out.emp_M = out.part_M/out.part;
    out.emp_R = out.part_R/out.part;
    out.emp_C = out.part_C/out.part;
    
    out.average_wage_M = integral(@(w) w.*h_M(w),w_lb,w_ub)*1/out.emp_M*1/out.part;
    out.average_wage_R = integral(@(w) w.*h_R(w),w_lb,w_ub)*1/out.emp_R*1/out.part;
    out.average_wage_C = integral(@(w) w.*h_C(w),w_lb,w_ub)*1/out.emp_C*1/out.part;
    
    out.total_cons_M = output.result.solution.phase.integral(9);
    out.total_cons_R = output.result.solution.phase.integral(10);
    out.total_cons_C = output.result.solution.phase.integral(11);
    
    % average income by occupation
    out.avg_inc_M = out.L_M*out.Y_M/out.emp_M;
    out.avg_inc_R = out.L_R*out.Y_R/out.emp_R;
    out.avg_inc_C = out.L_C*out.Y_C/out.emp_C;

    
    %% bound check info
    out.check_lb_L_M = check_lb(1);
    out.check_lb_L_R = check_lb(2);
    out.check_lb_L_C = check_lb(3);
    out.check_lb_K_B_M = check_lb(4);
    out.check_lb_K_B_R = check_lb(5);
    out.check_lb_K_B_C = check_lb(6);
    out.check_lb_K_E = check_lb(7);
    out.check_lb_K_S = check_lb(8);
    out.check_lb_b = check_lb(9); 

    out.check_ub_L_M = check_ub(1);
    out.check_ub_L_R = check_ub(2);
    out.check_ub_L_C = check_ub(3);
    out.check_ub_K_B_M = check_ub(4);
    out.check_ub_K_B_R = check_ub(5);
    out.check_ub_K_B_C = check_ub(6);
    out.check_ub_K_E = check_ub(7);
    out.check_ub_K_S = check_ub(8);
    out.check_ub_b = check_ub(9);    
    
    
    %% analyze ranges
    tab.min_w_transf = output.result.solution.phase.time(1);
    tab.max_w_transf = output.result.solution.phase.time(end);
    tab.f_1 = out.f(1);
    tab.f_end = out.f(end);
    tab.f_M_1 = out.f_M(1);
    tab.f_M_end = out.f_M(end);
    tab.f_R_1 = out.f_R(1);
    tab.f_R_end = out.f_R(end);
    tab.f_C_1 = out.f_C(1);
    tab.f_C_end = out.f_C(end);
    tab.min_l_scaled = min(output.result.interpsolution.phase.control(:,1));
    tab.max_l_scaled = max(output.result.interpsolution.phase.control(:,1));
    tab.min_V_scaled = min(output.result.interpsolution.phase.state);
    tab.max_V_scaled = max(output.result.interpsolution.phase.state);
    tab.fval_scaled = output.result.objective;
    tab.L_M_scaled = output.result.solution.parameter(1,1);
    tab.L_R_scaled = output.result.solution.parameter(1,2);
    tab.L_C_scaled = output.result.solution.parameter(1,3);
    tab.K_B_M_scaled = output.result.solution.parameter(1,4);
    tab.K_B_R_scaled = output.result.solution.parameter(1,5);
    tab.K_B_C_scaled = output.result.solution.parameter(1,6);
    tab.K_E_scaled = output.result.solution.parameter(1,7);
    tab.K_S_scaled = output.result.solution.parameter(1,8);
    tab.b_scaled = output.result.solution.parameter(1,9);
    
    struct2table(tab)

    %% save results
    vars = rmfield(out,{'wtilde','w','u','l','c','mtax','y','tax','avg_tax',...
                        'wtilde_noninterp','w_noninterp','u_noninterp','l_noninterp',...
                        'l_orig','u_orig','eta_orig','wtilde_orig','w_orig',...
                        'f_M','f_R','f_C','f', ...
                        'h_M','h_R','h_C','h','F_M','F_R','F_C',...
                        'F','H_M','H_R','H_C','H',...
                        'h_M_cond_part','h_R_cond_part','h_C_cond_part',...
                        'h_cond_part','H_M_cond_part','H_R_cond_part','H_C_cond_part','H_cond_part',...
                       'eta','phi_tilde_w','pi_w'});

    vectors = keepfield(out,{'wtilde','w','u','l','c','mtax','y','tax','avg_tax',...
                        'f_M','f_R','f_C','f', ...
                        'h_M','h_R','h_C','h','F_M','F_R','F_C',...
                        'F','H_M','H_R','H_C','H',...
                        'h_M_cond_part','h_R_cond_part','h_C_cond_part',...
                        'h_cond_part','H_M_cond_part','H_R_cond_part','H_C_cond_part','H_cond_part',...
                       'eta','phi_tilde_w','pi_w'});
    
    
    vectors.q_B = vars.q_B .* ones(size(vectors.w));
    vectors.q_E = vars.q_E .* ones(size(vectors.w));
    vectors.q_S = vars.q_S .* ones(size(vectors.w));
    vectors.prod_fun_eta_R = vars.prod_fun_eta_R .* ones(size(vectors.w));
    vectors.prod_fun_kappa_R = vars.prod_fun_eta_R .* ones(size(vectors.w));
    
    dir_out = '../../../data/generated/matlab/OptimTax/'
toc    
end

function out = welfare_term_foc(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:).*var.dF_dK_B(:));
end

function out = welfare_term_OB(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = fun_at_nodes * self.weights(:);

end

function out = welfare_term_B_1(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.f(:));
end

function out = welfare_term_B_2(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.dF_dK_B(:));
end

function out = welfare_shift_M(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_M(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:)./var.pi_w(:) .* var.Pi_M(:));
end

function out = welfare_shift_R(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_R(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:)./var.pi_w(:) .* var.Pi_R(:));
end

function out = welfare_shift_C(self, var)
% numerically integrate over participation costs
% integration is done of the transformed variable chi(phi),
% which transforms the domain of integration from
% [self.phi_lb, var.phi_tilde_w] to [-1,1] and transforms
% phi to chi accordingly. In contrast to welfare in GPOPSFunctions,
% no multiplication with var.wscale in the last row, since here I
% will integrate over w's, not over log(w)
    
    n_nodes = length(self.nodes);

    chi = @(phi,lb,ub) ...
          tools.ChangeOfVariables.change_variable_minone_one(phi,lb,ub);
    
    dchi_dphi = @(phi,lb,ub) ...
        tools.ChangeOfVariables.change_variable_minone_one_deriv(phi,lb,ub);            

    fun = @(u,phi,lb,ub) Psi_prime(self,u - chi(phi,lb,ub)).*self.part_dist.g_C(chi(phi,lb,ub)).* ...
          dchi_dphi(phi,lb,ub);
    
    nodes = self.nodes(:)';

    phi_tilde = var.phi_tilde_w;
    phi_tilde(phi_tilde < self.phi_lb) = self.phi_lb;
    
    u_mat = repmat(var.u(:),1,n_nodes);
    lb_mat = self.phi_lb.*ones(size(u_mat));
    ub_mat = repmat(phi_tilde(:),1,n_nodes);
    nodes_mat = repmat(nodes,length(phi_tilde(:)),1);

    fun_at_nodes = fun(u_mat,nodes_mat,lb_mat,ub_mat);

    out = ((fun_at_nodes * self.weights(:)).* ...
           var.V_prime(:)./var.pi_w(:) .* var.Pi_C(:));
end

function out = Psi_prime(self, V)
    out = V.^(-self.r);
end

function out = Psi(self, V)
    out = (V.^(1-self.r)-1)./(1-self.r);
end